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ABSTRACT 


Cullucation  schemes  are  presented  for  solving  linear  fourth  order  differential  equations 
in  one  and  two  dimensions.  The  var^iational  formulation  of  the  model  fourth  order  problem 
is  discretized  by  approximating  the  integrals  by  a  Gaussian  quadrature  rule  generalize<l  to 
include  the  values  of  the  derivative  of  the  integrand  at  the  boundary  points,  ('ollocaliun 
schemes  are  derived  which  are  equivalent  to  this  discrete  variational  problem.  An  efficient 
preconditionirr  based  on  a  tow  order  finite  differeoc  approximation  to  the  same  differential 
operator  is  presented  The  corresponding  multi-domain  problem  is  also  considered  and 
interface  conditions  are  derived.  Pseudospectral  approximations  which  are  C*  continuous  at 
the  interfaces  are  used  in  each  subdomain  to  approximate  the  solution.  The  approximations 
are  alscr  shown  to  be  f’’  continuous  at  the  interfaces  asymptotically.  A  complete  analysis 
of  the  collocation  scheme  for  the  multi  domain  problem  is  provided.  The  extension  of  the 
method  to  the  biharmonir  equation  in  two  dinwuisions  is  discussed  and  results  are  presented 
for  a  problem  defined  in  a  non  rectangular  domain 
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1  Introduction 


Spectral  metho<is  are  characterized  by  the  representation  of  the  solution  to  a  differential 
equation  in  terms  of  a  truncated  series  of  smooth  global  functions  which  are  known  as  trial  or 
basis  functions.  The  basis  functions  are  usually  chosen  to  be  the  eigenfunctions  of  a  singular 
Sturm  Liouville  problem  (Gottlieb  and  Orszag,  1977).  It  is  this  choice  which  is  responsible 
for  the  superior  approximation  properties  of  spectral  methods  over  other  standau’d  methods 
of  discretization.  For  linear  problems  possessing  smooth  solutions  these  eigenfunctions  yield 
expansions  that  converge  asymptotically  faster  than  any  finite  power  of  A'~'. 

Two  areas  of  research  in  spectral  methods  which  are  receiving  much  attention  at  the 
current  time  are  the  construction  and  analysis  of  well-posed  approximations  to  the  Stokes 
and  Navier-Stokes  equations  and  the  development  of  methods  which  can  be  applied  easily 
to  problems  defined  in  complex  domains.  With  respect  to  the  first,  it  is  well-known  that 
in  the  primitive  variable  formulation  the  velocity  and  pressure  approximation  spaces  need 
to  be  conipa^ible  to  avoid  problems  of  ill-conditioning.  This  is  similar  to  the  Ba^ulka- 
Brezzi  condition  required  for  the  corresponding  finite  element  approximation  spaces.  In  two 
dimensions  it  is  possible  to  avoid  this  difficulty  by  reformulating  the  governing  equations  in 
terms  of  a  stream  function.  The  governing  equation  is  then  fourth-order,  and  nonlinear  in 
the  case  of  the  Navier-Stokes  equations.  In  this  paper  we  seek  to  construct  pseudospeclral 
approximations  to  fourth  order  differential  equations  with  the  ultimate  goal  of  applying  them 
to  .solve  the  nonlinear  stream  function  formulatioo  of  the  Navier-Stokes  equaticms. 

Secondly,  the  development  of  techniques  for  handling  complex  geometries  is  essential 
if  spectral  methods  are  to  be  applied  to  problems  defined  in  mf>re  than  just  the  simplest 
domains.  The  basic  idea  behind  domain  decomposition  is  to  break  up  the  df>maiD  into  smaller 
simpler  subdomains  in  which  spectral  approximations  can  be  used.  The  approximations  are 
.suitably  linked  by  appropriate  interface  continuity  conditioas.  The  way  in  which  this  is 
implemented  is  important  if  the  full  power  of  the  spectral  method,  in  terms  of  the  accuracy 
of  the  approximation,  is  to  be  achieved. 

In  this  paper  we  shall  restrict  ourselves  to  the  model  fourth-order  problem  in  one  and 
two  dimensions.  Starting  from  a  variational  formulatioo  of  the  problem  we  shall  derive  a 
corresponding  collocation  problem  complete  with  interface  conditions.  In  a  domain  decom¬ 
position  setting  this  approximation  will  be  chosen  to  be  C*  continuous  implicitly.  In  addition 

continuity  across  the  subdomain  boundaries  is  achieved  asymptotically  as  the  order  of 
the  approximation  is  increased. 

Although  there  are  many  applications  of  spectral  methods  to  solve  second-order  elliptic 
partial  differential  equations  in  the  literature  there  is  tittle  previous  work  on  fourth-order 
problems  even  though  the  regularity  of  the  solution  to  these  problems  is  generally  higher  than 
for  second  order  problems.  Some  interesting  ideas  are  proposed  in  the  works  of  Moixboisne 
( 1984)  and  Orszag  ( 1971).  Bernardi  and  Maday  (1988)  give  a  survey  of  strategies  that  may 
he  employed  for  fourth-order  problems. 

Maday  and  Metivet  ( 1986)  have  studied  Chebyshev  spectral  and  pseudospeclral  approx¬ 
imations  of  the  stream  function  formulation  of  the  Navier-Stokes  equations.  They  prove  the 
convergence  of  the  schemes  and  derive  error  estimates  in  weighted  Sobolev  spares.  Kara- 
georghis  and  Phillips  ( I989a.l989b)  use  a  spectral  collocation  strategy  to  solve  for  the  laminar 
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flow  through  a  channel  contraction  again  using  a  stream  function  formulation  for  moderate 
values  of  the  Reynolds  number.  They  use  a  domain  decomposition  method  to  subdivide 
the  flow  region  into  rectangular  subdomauns  and  patching  to  piece  the  solutions  together,  in 
some  sense,  across  the  subdomain  interfaces. 

In  a  collocation  method  the  choice  of  the  coUoc^ion  points  is  crucial.  In  spectral  methods 
they  are  always  chosen  to  be  the  nodes  of  a  Gaussian  quadrature  rule  principally  for  two 
reasons.  First,  the  Lagrange  interpolating  polynomial  which  interpolates  data  at  these  nodes 
has  good  approximation  propertioi.  Secondly,  the  collocation  method  may  be  shown  to  be 
equivalent  to  a  variational  fcmnulation  of  the  problem  when  the  same  Gaussian  quadrature 
rule  is  used  to  approximate  the  integrals  appearing  in  this  formulation.  For  second-order 
problems  the  Gauss- Lobatto  nodes  are  used  because  the  boundary  conditions  can  then  be 
imposed  efficiently.  This  leads  to  an  optimal  error  in  the  resulting  spectral  approximation 
(Canuto  et  al.  (1987)).  For  fourth-order  problems  two  boundary  conditions  are  imposed  on 
the  solution.  1  hese  are  usually  of  Dirichlet  and  Neumann  type.  The  imposition  of  these 
boundary  conditions  is  facilitated  by  the  construction  of  a  generalized  Lagrange  interpolating 
polynomial  which  interpolates  the  function  at  the  interior  nodes  and  the  function  and  its 
derivative  at  the  boundary  nodes.  The  generalized  Gaussian  quadrature  rule  associated  with 
this  interpolating  polynomial  can  then  be  derived.  Quadrature  rules  of  this  form  are  quite 
well-known  in  the  theory  of  numerical  integratioo  (see,  for  example,  Golub  and  Kautsky 
(19^),  and  the  references  therein).  Golub  and  Kautsky  (1983)  describe  bow  the  weights 
in  these  quadrature  rules  may  be  determined  computationally.  In  this  paper  closed  form 
expressions  fur  the  weights  are  derived  using  the  properties  of  orthogonal  polynomials. 

We  show  that,  for  fourth-order  problems,  the  natural  choice  of  nodes  are  the  zeros  of  cer¬ 
tain  Gegenbauer  (or  uitraspherical)  polynomials.  Explicit  representations  for  the  quadrature 
weights  are  derived  for  evaluating  integrals  of  the  form 

J  ^  u>*(x)/(x)dr  . 

where  the  weight  function  takes  the  form 

Wx{i)  =  (I  -  X*)*  ,  A  >  -I. 

The  particular  form  of  these  weights  is  given  when  A  =  0  (the  Legendre  weight  function)  and 
A  =s  ~|/2  (the  Chebyshev  weight  function).  The  interior  nodes  in  the  case  when  A  =  -1/2 
are  the  wros  of  7^(x)  whereas  the  interior  Gauss- Chebyshev- Lobatto  nodes  are  the  zeros  of 

7"v.,(*). 

A  collocatioo  scheme  for  solving  a  fourth-order  model  problem  is  derived  by  omsideriog 
a  variational  formulation  ol  the  boundary  value  problem  with  suitably  denned  inner  prod¬ 
ucts.  The  two  formulations  are  shown  to  be  equivalent  if  the  inner  product  in  the  discrete 
variational  problem  is  defiMd  by  the  generalized  Gauss  quadrature  rule.  The  linear  system 
of  equatkms  which  derives  from  this  coHocatum  scheme  is  ill-cooditioned.  The  condition 
number  of  the  coefficient  matrix  scales  like  0(A'*)  where  S  is  the  order  of  the  approxi¬ 
mation.  An  efficient  precemditioneT  for  this  system  based  on  a  low  order  finite  difference 
approximatioo  to  the  same  differential  operator  is  presented.  The  combination  of  general¬ 
ized  Gaussian  quadrature  rules  with  spectral  methods  has  also  been  proposed  by  Bernard! 
et  al.  (1990).  This  idea  is  extended  to  multi-domain  problems  in  the  present  paper.  Pseu- 
dospectral  approximatioos  which  are  C*  continuous  at  the  subdomain  interfaces  are  used  to 
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approximate  the  solution  in  each  subdomain.  The  discrete  variational  problem  enables  us 
to  derive  interface  continuity  conditions  which,  in  the  asymptotic  limit,  result  in  con¬ 
tinuous  approximations.  The  variational  formulation  is  used  to  provide  an  analysis  of  the 
collocation  scheme  for  domain  decomposition.  The  analysis  shows  that  the  pseudospectral 
approximation  is  optimal  in  the  sense  that  it  is  of  the  same  order  as  the  corresponding  error 
in  the  best  approximation.  Numerical  results  are  presented  in  which  the  usual  exponential 
convergence  behaviour  of  spectral  approximations  is  exhibited.  Finally,  the  extension  of 
the  method  to  two  dimensions  is  described  and  numerical  results  presented  for  a  number  of 
model  problems.  An  application  of  the  method  to  the  solution  of  the  biharmonic  equation 
in  a  non- rectangular  domain,  an  L-shaped  region,  for  which  standard  spectral  methods  are 
not  applicable  is  presented. 

2  Variational  Formulation  of  the  Model  Problem 


In  this  section  we  consider  the  variational  formulation  of  the  fourth-order  model  problem. 
Consider  the  fourth-order  boundary-vadue  problem 

<t*u 


u(±l)«0. 


du 

dx 


(±1)  =  0, 


(1) 


where  f{x)  is  a  given  source  function.  It  is  well-known  that,  for  any  /  €  //~*(-I,I),  (1) 
has  a  unique  solution  u  €  1)  (  see  Grisvard  (1985),  for  example  ).  A  collocation 

scheme  for  solving  (1)  is  derived  by  considering  a  variatiooaJ  formulation  of  the  problem 
with  suitably  defined  inner  products. 

To  set  up  the  variational  formulation  we  need  to  define  function  spaces  for  each  A  >  -1. 
Let  Ll(-\,  I)  be  the  Hilbert  space  defined  by 


/»/  I  R  »  nreeasurable  ;  1 

=  \  /‘,u,*(xMx)dx<oo  / 


endowed  with  the  inner  product 


(u.»)a*  j  ^wx{x)u{x)v{x)dx. 


Vtft  also  introduce  the  Sobolev  space  1)  defined  by 


with  corresponding  fmrm 


» » ( x:  )■'■ 


(2) 
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Let  1)  be  the  subepace  1)  defised  by 

HU-U)^{v€Hli-Ul):  vi±l)  =  0;  «;'(±1)  =  0  }. 

Cmisider  new  the  btUnear  Imto  aA(-,>)  defined  on 

by  ^ 

OA(tt,»)  =  J  ^u"ix)[v>xix)v{x)Ydx.  (3) 

For  any  /  €  //~’(— 1, 1)  the  fourth-<vder  model  problem  (1)  is  equivalent  to  the  variational 
probknn  :  Find  u  €  1, 1)  such  that 

axiu,v)  =  (f,v)x  ,  ^v€  Hljai-lyl).  (4) 

Beroardi  and  Maday(I99l)  have  proved  the  following  result: 

Proposition  2.1  Let  X  satisfy  -1  <  A  <  1.  The  bilinear  form  o*  is  elliptic  on  W*_o(~l*  0 
and  eontinuoas  on  1)  x  1),  i  e. 

“  8“  iJjiS  «»(“.“) .  (s) 

Vu€  «;(-!, I),  .eHlj,(-i.i),  (6) 

respectively,  where  a  and  0  are  positive  constants. 

This  is  a  generalisation  ed  an  earlier  result  of  Maday  (1990)  for  A  s  -1/2.  An  immediate 
consequence  of  Propositioo  3.1  is  the  following  theorem: 

Theorem  2.1  Let  A  satisfy  -I  <  A  <  1.  For  any  f  €  Hx^{-1,\)  the  vari^ionai  problem 
(4)  has  a  snifse  sohtion  u  €  1).  Moreover,  it  satisfies 

II « Il3u<  C  B  /  .  (7) 


3  PMudcMp«ctral  AppraximalioD 

We  consider  the  pseodospectral  discretiaatkm  of  tbe  fourth-order  problem  (1).  Let  1) 

denote  tbe  space  of  algebraic  polynomials  of  d^rce  N  or  less  on  the  interval  (-1,1).  Let 

distinct  points  in  tbe  interval  (-1,1)  with  X|  =  -1  and  x^_i  =  1. 
Throughout  this  paper  we  take  /V  >  4  in  order  to  have  at  least  «ie  point  in  tbe  interior  of 
the  interval  (-1.1).  Suppose  that  tlw  values  of  some  function  /(x)  are  given  at  tbe  points 
.  1  <  i  <  ^  -  1.  together  with  the  values  /J  and  fff^,,  of  f'(x)  at  x  =  X],  and  x  =  x/r-t, 
respectively.  To  set  up  tbe  pseodospectral  approximation  of  ( 1 )  which  automatically  satisfies 
the  boundary  conditioos  it  is  necessary  to  coostruct  the  Lagrange  polynomials  for  this  data. 
Define  the  p<^ynomials  o(x)  and  IT(x  )by 

.(i)  =  (I  -  m*)  =  n%  -  X,)  (8) 

/•a 
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It  can  be  verified  that  the  La^ange  polynomial  of  degree  N  which  interpolates  this  data  is 
given  by 

Pn{x)  =  S  fA{x)  +  fiKiix)  +  fs-i^N-iix),  (9) 


where 


kj{x)  = 


_ nW*) 


2<j<N-2, 


T  _  "(*)  Si(») 

-  ftw  ■ 


The  corresponding  integration  rule  based  on  these  points 

N-l 

/  u?A(x)/(x)dx  =  53  tojfixj)  +  Wif{xi)  +  Ws-ifiiN-i)  , 

j-i 


can  be  shown  to  be  exact  for  all  /  €  Pts-ii-l,  1)  if  the  interior  nodes  Xj,2  <  j  <  N  -  2 
are  chosen  to  be  the  zeros  of  the  Gegenbauer  or  ultraspherical  polynomial  ^f  degree 

N  -  3.  The  location  of  the  nodes  are  determined  by  Newton’s  method  and  polynomial 
deflation.  For  the  sake  of  generality  we  consider  the  general  case  A  >  -1  here  although 
when  we  investigate  the  solution  of  differential  equations  we  will  only  consider  the  case  A  =  0. 
important  properties  of  the  ultraspherical  polynomials  G^*'*(x)  are  given  in  the  Appendix 
(hereinafter  the  reference  (A.m)  will  be  used  to  denote  equation  (m)  from  the  Appendix 
(rn  s  1,2, . . .)).  The  weights  depend  on  A  and  this  will  be  assumed  in  the  following. 

The  polynomials  ^^(x),  I  <  ;  <  iV  -  1  and  ^,(x),  ji  =  1,/V  -  1  defined  by  (10)  and 
(II),  respectively,  form  a  basis  for  P/v(“l.  !)•  Therefore,  choosing  /(x)  to  be  each  of  these 
polynomials  in  turn  we  obtain  explicit  expressions  for  the  S  +  I  weights; 


Wj  =  J  ^  Wx{x)hj{x)dx  ,  \  <  J  <  N  -  \, 

(14) 

Wf  =  j  ^wx{x)%,(x)dx  ,  j  =  l,S-\. 

(15) 

Although  only  the  boundary  weights  are  of  relevance  to  the  present  paper  in  that  they 
are  required  for  the  statement  of  the  multi-domain  collocation  problem  we  ^ve  details  here 
of  the  representations  for  the  interior  weights  as  well.  These  are  necessary  if  one  was  to  solve 
the  discrete  variational  problem  without  restating  it  as  a  collocation  one.  The  advantage  of 
doing  this  is  that  it  results  in  a  symmetric  system  of  linear  equations  to  be  solved  for  the 
unknown  nodal  values  of  the  solution. 

We  are  able  to  derive  an  original  result  in  which  explicit  representations  for  the  weights 
(9, 10)  are  obtained  using  the  properties  of  the  ultraspherical  polynomials  (see  the  Appendix). 
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Let  us  begio  with  the  weights  Wf,2  <  j  <  iV  —  2,  associated  with  the  interior  points.  The 
pdynomiab  TT(x)  and  are  related  by 

n(.)  =  ^.  (16) 

since  they  are  of  the  same  degree  and  have  the  same  zeros,  where  An-s  is  the  leading 
coefficknt  of  Thus  using  (10)  and  (14)  we  may  write 


Wj  = 


i: 


^>a(x)u(x)G}^!|>(x) 
(»  ~  *j) 


dx,  2<j<N  -2. 


(17) 


In  order  to  determine  the  value  ot  this  inte^al,  we  make  use  of  the  Christoffel  Darboux 
identity: 


Ik 


^N-i(x  -  y) 


An- 


(18) 


wiMre  'ifh,  (1  <  k  <  yV  -  3)  is  defined  by  (A.4).  Now  we  replace  y  by  x,  where  Gs-l\*))  ==  ^ 
then  (18)  redtKes  to 


7w-sAw.a  (» -  *,)  is 


If  we  now  multiply  both  sidei  of  (19)  by  (1  -  x*}^^’G0^*^(x)  and  integrate  with  respect  to 
X  over  (-1.1)  then  using  the  orthogonality  propvty  (A.3)  we  obtain 


-f-»  (*  - »,) 


which  enables  us  to  write 


2  <  >  <  yv  -  2, 


m 


since  g4^’^(x)  is  a  constant.  Unng  the  recurrence  rdatioo  (A.6)  with  n  s>  yV  -  3  we  write 
(17)  in  the  form 


2»^*r(yv>A-i)r(yv4.A) _ i _ 

'  (AT  - 3)!r(/V  +  W  +  2)  (|-,«)>G<;«''(x,)Gi;!.V)  ' 

for2<></V-2- 

Representatimis  for  the  boundary  weights  1st  aud  ^h-%  found  using  the 

integrals  (A. 8)  and  (A. 9).  Ccmsider  S»|  which,  in  view  of  (11).  (12)  (15)  and  (16),  may  be 
written  in  the  form 


®i  = 


(22) 
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Now  5t(x)  =  (I  -  x)*(I  +  *)  and  ther^ore  5J(— 1)  =  4.  The  cooditioo  (A. 5)  enables  us  to 
write  (22)  in  the  form 


u;»  = 


(-l)^(^-3)!r(A-f3)  /» 


4r(iV-fA)  J-,'*  '  '  '  ' 

The  integral  in  (23)  may  be  evaluated  analytically  using  (A. 8)  to  give 

2"''^*r(A  2)r(A  +  3)(iV  -  3)! 


Wi  = 


Similarly  we  can  show  that 

2^''^'^nA  f  l)r(A  +3)(A'  -3)! 


(A  +  3)niV-|.2A  +  2) 


ns  +  2A  +  2) 


{(A  +  2)A''  +  (A  +  2)(2A-  l).V-(4A*  +  9A  +  3)}.  (25) 


and  also 


u/,v-i  =  u.'i  ,  =  -n'l 


in  the  special  case  A  =  0  the  (legenbauer  polynomials  coincide  with  the  Legendre  poly¬ 
nomials  since  u'o(x)  =  I  Lhe  boundary  weights  are  given  by 

8r>.\^  -  2.V  -  3) 

=  irv-i)(\-iT:^Frii  ■  *■' ' 

and  the  interior  weights  by 

"  (.V  .rj^v  +  Di.v  +2)*(|  „  X*jir4*’, (/,)!*  ‘ 

for  2  <  j  <  (V  -  2  This  form  for  the  interior  weights  i*  derived  using  (A  6)  and  (A  10). 

When  A  =  -  I /2,  the  (<egmbaurr  polynomials  are  multiples  of  the  f'hebyshev  polynomials 
Tn{x)  *  roi»(n  ro<i“*r)  In  this  case  the  boundary  weights  are  given  by 


W|  ss  WiV-l 


3»(3.V*  -  6A'  +  I ) 
lOfiV  -2)(;V  -  D.V  ‘ 


!9|  at  =* 


4(iV  -2)(.V  1),V 


and  the  interior  weights  by 


t(S  ~  DfiV  -  ly 


(1  -x,*)(r;.,(x,)j 


-  .  2  <  ;  <  .V  -  2 


Having  written  down  an  expression  for  a  generalised  pseudospectral  approximation  (9) 
and  determined  the  weights  in  the  corresponding  quadrature  rule  we  are  now  in  a  position 


to  write  dowD  tlie  discrete  problem  correspoodiog  to  (4).  The  discrete  variational  problem 
OOTespoodiag  to  (4)  is  :  Find  us  €  Pjv(— 1, 1)  H  1, 1)  such  that 

ax(ujv,v/v)  =  ,  V  v;*  €  Pn{-1,  1)  n  1),  (33) 

wime  the  bilinear  form  (.,  .)x^  is  defined  by 

{f.9h.d=  £  «'7/(*j)s(*i)4'u;i({/fl){i,)  +  (/i>)(x^_i))  +  ufi((/p)'(x,)-(/5)'(xA,_,)],  (34) 
and  Xj,  2  <  j  <  iV  -  2,  are  the  zeros  of  Gjv^s^x). 

Theorem  S.l  The  vanattonal  prohUm  (33)  t»  efaisaienf  to  the  fotlomnp  collocation  prob¬ 
lem:  Find  UN  €  Pjv{~l,  1)  n  1)  sarfi  that 

<*'*’(»,)  =  /(»,)  .  2  <  }  <  ^  -  2,  (35) 

wAcrv  *j.2<)<5/-2  art  Uu  tent  tf 

Proof.  The  left  hand  side,  axiuN.VN)  »  integrated  by  parts  twice  to  give  the  following 
equivalent  problem  ;  Find  un  €  Pjv(  - 1 » I )  H  Hljti  ~  1 . 1 )  such  that 

»  {f,VN)xA  .  (36) 

foraUv)v€  PW(-1. 1)  H 1) 

Now  a  basis  for  the  space  Piv(  - 1 , 1 )  n  HljX  - 1, 1 )  are  the  polynomials  h,{x)  ,2<j< 

iV  -2,  defined  in  ( 10).  These  are  used  as  test  functiocMi  in  (36).  Since  €  Pis-ei  - 1 . 1  )n 

^  •  U  and  the  quadrature  rule  ( 13)  is  exact  for  all  polynomials  in  Pts-A  - 1 , 1 )  we  have 

(uS^’.h,)*  *  w,u^N^{x,)  .  2<  ]  <  S  -  2.  (37) 

Further,  from  the  definition  (34), 

(/.A,)m- w,/(^,)  .  2<;<  iV-2.  (38) 

and  therefore  since  w,  >  0,  for2<^<iV“2.  we  obtain  (35)  whkh  completes  the  proof 
of  the  theorem. O 

We  have  an  analogous  result  to  Theorem  2. 1  for  the  discrete  problem: 

Theorem  3.3  Let  \  »ati»fy  -  I  <  A  <  I  For  amp  /  €  r’®(-l,  I).  the  problem  (35)  ha$  a 
satfur  Mtfiiiteis  us  €  P/v(-l,  1)  n  1) 

Brrnardi  and  Maday  (lS19t)  establish  the  followiag  error  estimate: 

Theorem  3.3  Let  X  eahefp  -  I  <  A  <  1  ff  the  eolmtion  u  of  (39)  brlonp*  to  for 

a  not  number  er  >  I,  sad  if  the  data  f  is  imch  that  the  function  (I  -  x*)*/  belongs  to  a  »pact 
/fj(-l,  I)  for  a  real  number  p  >  1/2.  the  foUoennp  error  rjiitmaif  between  the  .*olutwn»  of 
(29)  and  (33)  m  $att»^ed 

\\u-UN  h.y<  CiS*-*  [i  u  iU.»  (I  (I  -  x*)t/  11,.,). 

n 


(39) 


The  collocation  method  (38)  results  in  a  system  of  equations  for  the  values  uj  of  ua^(i) 
at  the  points  Xj  ,2  <  j  <  N  —  2.  The  pseudospectral  collocation  approximation  is  then 
given  by 

N-i 

=  XI  >  (40) 

J=t 

where  hj(x)  is  given  by  (10)  (cf.  (9)). 


The  generalization  of  the  collocation  method  (35)  for  problems  with  inhomogeneous 
boundary  conditions  is  straightforward.  The  nature  of  the  pseudospectral  approximation 
(9)  is  such  that  inhomogeneous  boundary  conditions  are  satisfied  exawrtly  by  simply  insert¬ 
ing  the  specified  values  directly  into  (9).  If  1)  is  the  subspace  of  1)  which 

consists  of  those  functions  that  satisfy  the  given  inhomogeneous  boundary  conditions  then 
we  have  the  collocation  problem  :  Find  us  6  /V(-l,  1)  H  //*  g(-l,  1)  such  that 

“&'’(*j)  =  /(*j)  ,  2  <  j  <  yV  -  2, 

where  x,  ,2  <  j  <  N  —  2  are  the  zeros  of 

4  Preconditioning 

The  collocation  problem  (35)  can  be  restated  in  the  form  of  a  linear  system  of  algebraic 
equations 

=  b  (41) 

where  u  is  the  vector  of  values  of  u/v(x)  at  the  collocation  points  X;,  2  <  ;  <  N  —  2,  b  is  the 
vector  of  values  of  /(x)  at  these  points  and  A  is  the  {N  —  Z)  x  (N  —  3)  matrix  whose  entries 
are  defined  by 

a1"'’(x,).  2<j,k<N-2. 

The  fourth-order  pseudospectral  differentiation  operator  A  has  postive,  real  eigenvalues.  The 
extreme  eigenvalues  of  A  are  shown  in  Table  1.  In  this  table  we  sec  that  the  largest  eigenvalue 
of  A  .scales  like  N*  while  the  smallest  eigenvalue  is  independent  of  /V.  Therefore,  since  the 
condition  number  of  <4  is  0{S*)  an  efficient  preconditioner  is  necessary  for  the  accurate 
inversion  of  (41). 

Orszag  ( 1980)  proposed  a  preconditioner  for  spectral  methods  based  on  a  low-order  finite 
difference  approximation  to  the  same  differential  operator.  The  advantages  of  such  a  pre¬ 
conditioner  are  that  it  is  sparse,  easily  invertible  and  yields  an  inverse  close  to  the  inverse 
of  the  original  spectral  operator.  Therefore  we  propose  using  a  second-order  finite  difference 
operator  as  our  preconditioner.  This  requires  the  solution  of  a  pentadiagonal  system  which 
may  be  performed  very  *'Hiciently  and  stably  using  Gaussian  elimination. 

For  3  <  i  <  iV  —  3  the  second-order  finite  difference  approximation  to  u(x)  at  the  point 
X,  is 

u<"'’(x,)  sa  o,u(x,.j)  6,u(x._,)  c,u(x.)  -f  d,u(x,+i)  +  e,u(x,+2) 

where 

24 

<1^  22  T,,,  I  ■  ,1  . . .  . 

(x,.a  -  x,){x,_a  -  x._i)(x,_a  -  x,+,)(x,_3  -  x.+a)’ 
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_ -24 _ 

(Xi— 1  Xi)(Zt-.2  Zi_i)(x,'_i  Xt4.i)(Xi-i  ^t+2) 

_ -24 _ 

(Xj  3  35j^l)(Xi_i  Xi^.i)(Xi^l  ^t+s) 

_ 24 _ 

(Xj  Xi+2)(^«— 3  3Cj4.2)(x,'_j  Xt^2)(^(+1  ®«+2) 

Ci  =  — (oj  +  6i  +  <i, ■  +  c,). 

A  similar  formula  holds  at  X2  and  xs-2  after  taking  into  account  the  homogeneous  Neumann 
boundary  conditions.  Let  H  denote  the  finite  difference  differentiation  matrix  defined  by 
the  above  equations.  We  are  interested  in  the  eigenvalue  spectrum  of  the  operator  H~^A 
since  this  governs  the  rate  of  convergence  of  the  preconditioned  iterative  method  for  solving 
(41).  The  eigenvalues  of  H~^A  are  real  and  positive.  The  extreme  eigenvalues  of  H~^A 
are  shown  in  Table  2.  Again  the  smallest  eigenvalue  remains  independent  of  the  choice 
of  N  while  the  largest  eigenvzJue  grows  very  slowly  with  N.  The  entries  in  this  table 
demonstrate  the  effectiveness  of  //  as  a  preconditioner  for  A.  Haldenwang  et  al  (1984) 
showed  theoretically  that  the  eigenvalues  of  the  corresponding  preconditioned  second-order 
pseudospectral  differentiation  operator  lie  between  1  and  (7r/2)^.  From  this  result  one  would 
expect  that  the  eigenvalues  in  the  case  of  the  fourth-order  problem  to  lie  between  1  and 
()r/2)^.  We  can  see  from  Table  2  that  they  do  indeed  lie  between  these  bounds. 

5  Analysis  of  the  Multi-Domain  Problem 

Given  a  fixed  integer  Af  we  consider  a  partition  of  (—1, 1)  into  M  subintervals  /„  where 

An  —  (dn»>^tn+l)> 

and  the  dm  are  M  A  \  points  in  (—1, 1)  such  that 

-1  =  do  <  d\  <  •  •  •  <  dut-i  <  dM  =  I- 

Associated  with  each  subinterval  /„,  we  define  a  set  points  xJ*  ,  1  <  j  <  N  and  weights 
.  1  <  j  <  ^-1  =  1,JV-1,  which  correspond  to  a  generalized  Gaussian  quadrature 

rule  of  the  form  (13)  defined  on  /„.  Let  h^{x)  1  <  j  <  N  —  I,  and  V^{x),  j  =  l,N  —  l,he 
the  corresponding  interpolating  functions  which  have  compact  support  on  the  interval  Im- 
We  introduce  the  finite  dimensional  spaces 

F/v  =  €  L’(-l,  1) :  (j>\i„e  FW(/m)}, 

where  N  is  some  integer  and  FW(A)  denotes  the  set  of  all  polynomials  of  degree  less  than  or 
equal  to  N  over  A.  In  order  to  discretize  the  space  //q(— 1, 1),  let  us  introduce  the  spaces 

XN  =  YNnHU-lA\ 

L’(-i,  1) :  <i>  ke  PNiU)  n 

The  elements  of  X/v  are  continuous  and  have  continuous  derivatives  at  the  points  dm,  1  < 
m  <  M  —  I,  and  vanish  along  with  their  first  derivatives  at  x  =  ±1. 


bi  = 
di  = 
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In  this  paper  we  shall  only  consider  the  case  A  =  0  as  far  as  domain  decomposition  is 
concerned.  This  is  the  only  value  of  A  for  which  the  weight  function  over  (  —  1,1)  is  the 
same  as  the  weight  function  over  each  of  the  subintervals  /„,  \  <  rn  <  M .  Throughout 
this  section,  in  which  A  =  0,  the  zero  subscript  has  been  deleted.  For  example,  a(.,.)  is 
used  synonymously  with  a(.,.)o  and  the  corresponding  norm  notation  has  also  been  altered 
accordingly.  We  now  define  the  discrete  problem:  Find  ujv  fc  such  that 

a(uN,v.v)  =  (/, yyv)M,  V  Vs  €  Xs,  (42) 

where  the  bilinear  form  (.,  .)m  is  defined  by 

(/.»)«=£(/•«)-  («) 

m  =  l 

where 

U-gU  =  z  +  <"ri(/9)(^r)  +  (/9)(^5-,)I  +  KK/gruT)  - 

Lemma  5.1  For  any  rtal  number  <r  >  2  and  for  any  o  €  Hq(  —  \.  1)  H  I)  we  have 

II  ♦  -  ll>  <  CiV'-'  II  «  II,. 

where  is  the  orthogonal  projection  operator  from  Hq{-1,  1)  onto  Fs(-\,  \)C\Hq(-\.  1). 
Proof.  See  Bernardi  and  Maday  (1991).  □ 

Since  //*(-!,  1)  is  contained  in  C‘((-I,  I))  we  can  show  that  for  any  €  //*(-!,  1).  there 
exists  a  cubic  polynomial  ^  such  that  ^  ^  €  Ho[~\,  1)  and  for  any  real  number  .s  >  0. 

II  Oo  IU<  C  II  ^  Ha  . 

Now  define  an  operator  wjf  by  -  ^)  +  ^  from  //^(-l.  1 )  onto  Psi-\,  !)■  So 

that  if  ^  €  //'(  —  1, 1)  then  by  L^ma  4.1, 

II  ^  iii  =  II  Ha 

<  c/v’-' IM  -  A)  lU 

s  C/V’-"|UIU- 

We  can  easily  verify  that  this  operator  satisfies 

(x^^)(±l)  =  ^d:l),(x»,^)'(±l)  =  ^'(il). 

Theorem  5.1  There  exials  an  operator  iff  from  Hq{~\,  1)  onto  Xs  satisfying 

ll»^-*^0||a<CN*-'IN'IU.  (44) 

for  any  function  €  //'(—1, 1)  n  1)  with  c  >2. 
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Proof.  We  recall  that  for  a  general  interval  (a,  6)  there  exists  a  projection  operator  icn  from 
b)  onto  Piv(o,  6)  satisfying 

II  ti;  -  |Ua(,.>,<  C II  u;  !!>/.(,»,  (45) 

7rfi/w(a)  =  ui(a),  (»^u;)'(o)  =  it;'{a),  (46) 

vsw(b)  =  u;(6),  (»nu>)'(6)  =  w'{b),  (47) 

for  all  uj  €  H‘'{a,b). 

Let  us  define  the  projection  operators  r/v.m.  for  1  <  ni  <  Af,  as  being  the  projection 
operators  from  onto  P\(/»).  We  deduce  that  the  element  defined  on  each  /„ 

by 

»N^(-f)  =  *'iV,m0(-r),  V  X  €  /m, 
is  an  element  of  fV/( - 1 , 1 )  H  Hq(-  1,1)  that  satisfies  due  to  (48) 

II  ^  l|j<  II  0  lU  .□  (48) 

Define  J/v/  to  be  the  Lagrange  interpolating  polynomial  which  interpolates  the  function 
/  at  the  iV  —  3  interior  collocation  points  of  the  generalized  Gauss  quadrature  rule  on  (  —  1,1). 
Then  Bernardi  and  Maday  (1991)  have  shown  that 


Lemma  5.3  For  any  rtal  number  p  >  1/2  and  for  any  0  such  that  the  function  (1  —  x^)t0  € 
1),  Me  following  tneguoHty  holds 


II  (I  -  -  JnI)  ||o<  CN''’-'  il  (I  -  i’)t/  II.  .  (49) 

Lemma  5.3  For  any  rtal  number  >  5/2  and  for  any  f  such  that  the  function  |(dm+i  “ 
s)(i-d^)\yy€H'-{U),  then 


sup 


(/.  '’s)u  ~  if,Vs)m 


II  VS  l|#f>(;«) 

where  (..  )/„  w  the  L*  inner  product  on 


<  CN‘'’-’~  II  1(4..,  -  i)(i  -  4,)|"'V  IU~„.,, 

(50) 


Proof.  The  generalized  Gauss  quadrature  rule  on  1^  »*  exact  for  any  polynomial  in 
/j/v-3(/m)  and  so  for  any  vjv  €  Ps{lw,)  H  we  have 


{f,Vs)u~if’Vs)m  ~  {f  -  Jsf,Vs)u, 


where  Js  is  the  Lagrange  interpolation  operator  at  the  N  -  .1  interior  nodes  of  a  generalized 
Gauss  rule  on  the  interval  l^-  We  recall  that  the  mapping  w  *-*  a’/((dm+i  —  x)(x  —  d„)]^  is 
continuous  from  //^(/m)  into  £,*(/»,).  Then  we  ran  write 


(/■  Vs)u  -  (/.  l'(v)m  <  G*  II  ((<C>I  -  •?)(•*■  -  <fm)j’'^*(/  -  Jsf)  IIl»(/„)II  Vs  • 


Finally  using  Lemma  4.2  we  obtain 

(/.i-n);, -(/.w),  <  II  1(4..,  -  i)(x  -  4.)|«V  IIh~„.i1I  tJN  |l«.(,.,,  (51) 

from  which  we  deduce  the  result.  □ 
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Theorem  5.2  Let  us  suppose  that  the  solution  u  to  (32)  belongs  to  H"(  —  \,\),  for  a  real 
number  <r  >  2  and  that  —  x){x  —  </«,))*/  €  //'’"*(/„)  for  a  real  number  Pm  >  1/2  for 

each  m  =  0,  —  1.  Then  the  following  error  estimate  holds: 

II  u  -  ||,<  C(/V>-’  II  u  II.  +  /V'"-'-  II  Ki..,  -  I)(I  -  (52) 

m=0 

Proof.  Let  us  define  u*  =  u  —  u®  and  u%  =  uyy  —  ti^,  where  u®  and  u%  are  piecewise  cubic 
polynomials  such  that  u*  |/„€  //o(/m)  and  u%  |/^G  PwC/m)  H  //o{/m),  0  <  m  <  Af  -  1. 
Then  Proposition  2.1,  together  with  (4)  and  (42),  gives  for  «iny  vs  €  Zn, 

l|n^-nAf||3  $  Cia{ul, -VN,u}j  ~vn) 

=  C{a{u’ -  vs,ul, -Vs)  -  {f,u’’fj -vs) {f,u'’fj -vn)m),  (53) 


from  which  we  obtain 

.1  .  •  >i  t  It  •  11  .  {f,^N)-{f,^s)M\ 

II  U  -^s  Il3<  'of  II  «  -  Ha  +  sup  - - - - - >  . 

{ys6Zs  y>NiZs  II  111  J 

We  choose  vs  =  and  use  Theorem  4.1  to  show  that 

inf  II  u-  -  u/v  H3<  CN^-^  II  u-  ||„  (T  >  2. 

Viv€^/v 

Since  ws  €  2s  we  have  on  each  interval  /« 

{/.«'/v)/«  -  (f,»fs)m  =  (/  -  JNf,Ws)u- 

We  can  also  show  that 


{f,ws)  -  (/,  ws)m  ^ 

II  .....  II-  - 


sup 


sup 


{f,^N)u  -  {fy^N)r 


||3  II  ||«>(/»,) 

and  therefore  using  Lemma  4.3  we  may  deduce  that 


s  c  E  II  ll«~«-i  • 


sup 

wft^Zs 


||3 


niaO 


(54) 

(55) 


(56) 


Since 

II  U  -  u^  ||3<1|  u*  -  u%  Ha  +  II  U®  -  lla 

and 

II  u®  -  u®^  ||3<  C  II  u-  -  u%  ||„ 

then 

II  «  -  u/v  !l2<  C*  II  -  u)v  II3  • 

Finally  using  (54)-(56)  we  obtain  the  result.  □ 

We  now  set  up  the  collocation  scheme  for  the  domain  decomposition  problem.  We  define 
Us  €  Xs  which  interpolates  data  at  the  points  x”  ,  1  <  j  <  yV  —  1,  1  <  m  <  Af  by 

«/v(x)  =  Z  u7A7(x)  +  (uW(x)  +  (u%_,K;_,{x)  ,  x€lm,  (57) 

j=i 
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where 


(58) 


ujj  =  ,  {u%  =  (u')r^‘  .  1  <  m  <  A/  -  1. 


Theorem  5.3  The  variational  problem  (4S)  with  the  discrete  inner  product  defined  by  (4S) 
is  equivalent  to  the  following  collocation  problem  :  Find  un  ^  Xs  such  that 

«N{^r)  =  /(*")•  2<j  <  A^-2.  1  <m<  A/,  (59) 

UN(*r'+)-«X(x;5_,-)  =  u,r*r(xr>+)  +  «,j)_,r(x;;;.,-), 

l<m<Af-l,  (60) 

=  -wr»r(xr>+)-ti,;5.,r(x;5-,-) 
-ror^v(xr‘+)  -  n5:?.,r'(x]5.,-), 
I<m<A#-l,  (61) 

and  where  r(x)  =  u)y(x)  —  /(x). 

Proof.  Let  lu  examine  a(ujv,  v/v)  defined  by  (3).  By  linearity  we  may  write  the  integral  on 
the  right-hand  side  of  (3)  as  the  sum  of  integrals  over  each  subinterval  /m  for  1  <  m  <  M. 
Subsequent  integration- by- parts  twice  gives 

a(uw,u/v)=  ^  /  u;;(x)vw(x)dx  -  {(«nVw](x;J.,)  -  (u>w)(x;5.,)},  (62) 

m«l  mm% 

where  \f](y)  s  /(y+)  -  /(y-)  denotes  the  jump  at  x  =  y  in  /. 

We  choose  as  our  basis  for  the  space  Xn  the  polynomials  A]^(x),  2  <  j  <  Af  —  2, 
1  <  m  <  M  —  \  and  Aj5_i(x),  ^_|(x),  1  <  m  <  Af  —  1.  The  use  of  these  polynomiads  as 
test  functions  in  (42)  with  the  discrete  inner  product  given  by  (43)  results  in  (59)-(61)  which 
completes  the  proof  of  the  theorem.  □ 

Remark  1  Note  that  in  view  of  the  expressions  for  the  wei^Us  given  in  (27)  and  (28), 
u;^  =  ws  =  0{N~^) ,  =  —TSu  =  0(N~*) ,  as  N  -*  oo  , 


and  therefore  from  (60)  and  (61)  we  can  write 

■-;(xr'+)-«;(*5-.-)  =  o(Af->), 

«;(»:*'+)- «sw-.-)  =  o(w-') , 

as  N  —*  oo.  Thus  we  have  second  and  third  order  continuity  at  the  interface  asymptotically, 
as  N  —*  oo. 


6  The  Biharmonic  Problem  in  Two  Dimensions 

Consider  the  biharmonic  problem 

VV(*,  y)  =  /(*.  y) ,  in  ft  ,  (63) 

V>(*,y)  =  yi(^.y) .  on  r  ,  (64) 

^(^.y)  =  y3(*.y) .  on  r ,  (65) 

where  ft  =  (—1,1)  x  (—1,1)  and  F  is  the  b6undary  of  ft.  Grisvard  (1985)  shows  that 
provided  the  boundary  data  satisfies  certain  compatibility  conditions  there  exists  €  //*(ft) 
satisfying  (64)  and  (65).  Since  we  are  primarily  concerned  with  the  domsun  d  ^position 
problem  we  only  consider  the  case  when  the  weight  function  is  unity.  The  ai  for  the 
single  domain  problem  is  thus  greatly  simplified. 

In  order  to  write  down  the  variational  formulation  of  the  problem  (63)*(65)  we  define  the 
bilinear  form  on  A/^(ft)  x  //^(ft): 


0(0,*)= yj(^(vv)(vv)<fa<i> . 


The  biharmonic  problem  (63)-(65)  is  then  equivalent  to  the  following  variational  problem; 
Find  Ip  €  //^(ft)  such  that  (0  —  0*)  €  Wo(ft) 


«(V»,0)  =  (/.0),  for  aU  0  € //o^ft) , 


where 


(/.^)  =  ' 


We  see  that  0  is  a  solution  of  the  variational  problem  (67)  if  and  only  if  0  =  0  —  0^  is  a 
solution  of  the  problem:  Find  0  €  Hq{U)  such  that 

a(0, 0)  =  (/,  0)  -  a(0*,  0) ,  for  all  0  €  .  (68) 

It  can  be  easily  verified  that  the  bilinear  form  o(.,.)  defined  by  (63)  is  continuous  and 
elliptic  on  Ho(n)  x  f/o(n)  and  hence  that  problem  (71)  has  a  unique  solution  in  Ho(ft)  for 
/  6  f/-’(ft). 

Let  /’/v(ft)  be  the  space  of  algebraic  polynomials  of  degree  at  most  N  in  each  co-ordinate 
direction.  The  collocation  problem  associated  with  (63)-(65)  is: 

Find  0/v  €  P/v(ft)  FI  //’(ft)  such  that 

v^0w(*.y)  = /(ic.y) .  (*.y)e//,  (69) 

0/v(*.y)  =  yi(*.y) »  (x,y)6  5uT,  (70) 

=  (r,»)€5ur,  (71) 

?^(*.ll)=  (*,»)€  r,  (72) 


where  dfdn  and  dfdt  represent  differentiation  pormal  and  tangential  to  F,  respectively,  the 
sets  A,  5  and  T  are  defined  by 


«  =  -  -2), 

S  =  {te,±l),(±l.fi):2<i<W-2}, 

T  =  {(±1,±1)). 

and  C»iJL3(6)  =  0,  2<i<yV  —  2.  There  are  a  total  of  (N  +  1)’  linear  equations  for  the 
(N  -f  1)^  unknowns.  The  dimension  of  Pn{0)  is  {N  1)’.  The  basis  functions  in  2D  are  the 
tensor  product  of  the  cme-dimensional  basis  functions  given  by  (10)  and  (11). 

We  define  the  two-dimensional  discrete  inner  product  in  an  analogous  way  to  (34)  by 
applying  the  quadrature  rule  in  each  co-ordinate  direction.  So  in  the  case  when  one  of  the 
functions  V'  or  ^  belongs  to  we  have 

=  EE  ■  («) 

im2  jm3 

Next  we  define  the  discrete  bilinear  form  a/v(.,.)  by 

aN(0.^)  =  (V^0,^)/v.  (74) 


Theorem  0.1  If  there  is  a  funeiion  satisfying  the  boundary  conditions 

(61)~(68)  then  the  collocation  problem  (69)-(72)  is  equivalent  to  the  variational  problem: 
Find  €  Aw(n)  n  //>(n)  such  that  (^/^  -  and 

0^(0//. ^n)  =  {f^4s)N  .  for  oil  ds  €  Pn{^)  H  /fo(n) .  (75) 

Proof.  On  each  horizontal  or  vertical  side  of  H,  0/v  end  0jv  are  polynomials  of  degree 
yV  satisfying  N  I  conditions  and  so  are  identical  on  F.  The  same  argument  applies  to 
their  normal  derivatives  and  so  (0^  -  0jv)  €  Pn{^)  H  //o(n).  If  we  now  choose  0N(x,y)  = 
kj{x)hit(y),  2  <  j,k  <  N  —  2,  then  (75)  implies  (69)-(72).  Conversely,  since  these  {N  —  3)^ 
polynomi^s  form  a  basis  for  Pn  n  //o(n),  (69)-(72)  implies  (75).  □ 

Let  us  now  turn  our  attention  to  the  problem  of  domain  decomposition,  and  for  simplicity 
restrict  ourselves  for  the  moment  to  the  case  when  0  is  divided  into  two  subdomdns  with 
interface 

7  =  {(0.y):  -i<y<i}. 

We  define 

n,  =  {(x,y):-l<x<0,  -l<y<l}, 
fta  =  {(x.y) :  0  <  z  <  1,  -1  <  y  <  1}  , 

and  F*  is  the  boundary  of  for  k  =  1,2.  Define  the  subspace  V  of  x  //^(Dj)  by 

V  =  {♦  =  (V-',0’)  €  H’(n,)  X  H’(n,) :  V-'  =  It’,  ^  ^  om)  , 
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and  the  subspace  Vq  of  V  by 

VJ,  =  {♦  =  €V  :  V'"  =  ^  =  0  on  r  for  m  =  1,2  }  . 

an 

We  consider  the  bilinear  form  defined  on  V  x  V  by 

.(♦,♦)  =  a, (0‘.*')  +  a,{^’,#'),  (76) 


where 

<>*(v!>‘,^‘)  =  /jf  _(VV*)(VV*)<fa</»  , 

We  can  show,  using  Green’s  theorem,  that  if  4  €  Vo  then  the  above  bilinear  form  may 
be  written  as 


(77) 


If  there  exists  4**  €  V  satisfying  (64)-(65)  then  the  variational  problem  is:  Find  4^  €  K  such 
that  4>  -  4**  €  Vo  And 


<’*•*)  =  f  L  L  .  (TO) 


where  /*  is  the  restriction  of  /  to  Os. 

The  variational  problem  (78)  is  equivalent  to  the  following  interface  problem: 


V**‘  =  /‘,  inU,,  t=1.2, 

(79) 

*‘  =  ji .  on  rnr» ,  t  =  1,2 , 

(80) 

—  =  onrnr*,  *  =  1,2, 

(81) 

,  ,a 

*  =*  '  8r  =  -8r’ 

(82) 

Define  the  finite  dimensional  space  Vn  by 

K/v  =  (♦  =  (0',^’)  €  P;v(ni)n  A/’(n,)  X  PAr(na)n  //’(Dj) : 


and  the  subspace  Vsjo  of  V/v  by 


dx 


on  7}  , 


Vn,o  =  {♦  =  €  V/v  :  0"*  =  ^  =  0  on  r  for  m  =  1,2  }  . 

In  the  case  when  one  of  4^  or  ^  belong  to  V/v,o  we  define  a  discrete  inner  product  by 
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where 


N-2N-7 

(!(’*. *‘)‘n  =  E  E 

iac3  js3 

i=3 

for  k  =  1,2,  where  d/dn  is  the  normal  derivative  to  ft*.  The  discrete  bilinear  form  on 
Vn  X  is  defined  by 

a/vC*, ♦)  =  a)y,(0*,  ^‘)  +  aji,(V>^,  , 

where 


for  ib  =  1,2. 


Theorem  6.2  If  there  is  an  element  €  Vn  satisfying  (80)  and  (81)  then  the  variational 
problem:  Find  '^n  €  V/v  such  that  (^n  —  ♦)v)  €  V/v,o  <*«</ 

««(♦«.♦«)  =  (/',#■)},  +  (/’.*’&  .  (83) 

for  all^  sz  €  Vsja,  ts  equivalent  to  the  collocation  problem:  Find  €  Vn  such  that 

VVjv(*,y)  = /*(*,y) ,  (z,y)eH*,  *=1,2,  (84) 

V'Mx.y)  =  yi(i,y) ,  (z,y)€5*ur*,  fc=l,2,  (85) 

=  y3(x,y) ,  (z,y)€5kU7k,  fc=l,2,  (86) 

?%(x.y)  =  ^(x.y) .  (i,y)€rfc,  fc  =  i,2,  (87) 

cmcrC  (ft 

-u,;i(vvi!,-/)(o,fi)i 

-  /)(o,f,)i 

-®;i^(VVi!,-/)(0,f,)),  2<j<N-2,  (88) 

-  V'JS<)(0,(y)  =  -  /)(0,{,)1 

+  «.»[(VV;,-/)(0,{y)l,  2<j<N-2,  (89) 
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when 


«*  =  {((i.ii)-i<i.j<N-2),k=l,2, 
S,  =  {({•.±l),(-l.(i):2<i<fl^-2), 
S,  =  {(f?,±l).(l.<i):2<i<Ar-2), 
r,  =  {(-1,±1).(0,±1)). 
r,  =  ((o,±i),(i.±i)}. 


fi=(fc-l)/2.  {.■=(l+6)/2.  2<i<N-2. 

Proof.  We  can  show  (♦«  —  ’I'w)  €  VVn  “  in  fhe  proof  of  Theorem  5.1 .  If  we  now  choose 
as  our  test  functions  the  following: 

+  l)Ai(y) ,  ^’(x,y)  =  0,  2<kJ<N-2, 

,  ^*(x,y)  =  0,  ^*(x,y)  =  M2j:- l)A|(y),  2<kJ<N-2,  .  . 

0‘(x,y)  =  hA,-,(2x  +  l)h,(y) ,  <f>\x,y)  =  h,(2x  -  l)h,(y) .  2<l<N-2, 

.  ^‘(x>y)  =  ^n-i(2x  +  l)h/(y) ,  ^’(x,y)  =  S,(2x- l)h/(y) ,  2<l  <  N  -2  , 

then  we  obtain  immediately  (84),  (85),  (88)  and  (89).  Conversely,  since  these  2{N -Z){N -2) 
test  functions  constitute  a  basis  for  Ysfit  (84)'(^)  implies  (83). 


7  Numerical  Results 

The  quadrature  rule  (8)  is  used  to  compute  approximations  to  the  integrals 

(a)  /i,  u;A(x)e*cois()ri)dx, 

(b)  /!,  ufA(x)xe*di, 

when  A  =  0  and  A  =  — 1/2.  The  errors  in  the  quadrature  rule  are  given  in  Tables  3  and  4  for 
integrals  (a)  and  (b),  respectively,  for  different'  values  of  N.  The  quadrature  rule  evaluates 
the  integrals  accurate  to  machine  precision  for  a  value  of  as  small  as  17. 

7.1  1-D  Problems 

Numerical  solutions  to  the  fourth-order  model  problem  (1)  are  obtained  when  the  exact 
solution  is  given  by 

(a)  u(x)  =  (1  -  x^)*sm(irx), 

(b)  u(x)  =  1  -f-  sm(2xx). 

In  example  (a)  the  boundary  conditions  are  homogeneous  whereas  for  (b)  we  have  inho- 
mc^eneous  boundary  conditions.  The  differential  equation  is  collocated  at  the  generalized 
Legendre  and  Chebyshev  nodes  given  by  the  zeros  of  (1  -  x*)f^_,(x)  and  (1  -  x’)rjv_,(x), 
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respectively.  The  error  in  the  numerical  solution  is  measured  using  weighted  norms  based 
on  the  corresponding  generalized  quadrature  rule.  The  infinity  norm  is  also  given  to  show 
the  maximum  pointwise  error  at  the  collocation  points.  These  are  displayed  in  Tables  5  and 
6  for  examples  (a)  and  (b),  respectively,  where  we  define 

II  «  h.^  =  IL  -  («/v-i)’)))‘^*. 


II  e  ||ao=  max  |  e.  |  , 

and 


ej  =  u{xj)-UN{ij),  j  =  -  \  , 

where  the  points  x,,  1  <  j  <  N  —  1,  are  the  generalized  nodes.  The  usual  exponential  con¬ 
vergence  of  spectral  approximations  to  smooth  solutions  of  differential  equations  is  observed 
with  accuracy  to  machine  precision  being  obtained  when  N  =  24. 

Next  we  apply  these  techniques  in  the  case  of  domain  decomposition.  For  simplicity  we 
consider  a  partition  of  the  interval  (-1, 1]  into  the  two  subintervals  [—1,0]  and  [0, 1]  with 
common  point  x  =  0.  We  solve  again  the  model  problems  (a)  and  (b)  using  the  collocation 
scheme  (59)-(61).  The  corresponding  error  norms  are  shown  in  Tables  7  and  8,  respec¬ 
tively.  The  mono-domain  and  two-domain  spectral  approximations  converge  exponentially 
as  expected.  The  two-domain  approximation  converges  slower  than  the  mono-domain  ap¬ 
proximation  for  the  same  total  number  of  collocation  points  since  for  the  problems  considered 
here  there  is  no  particular  advjtntage  to  be  gained  in  using  the  former  since  the  solutions 
are  smooth  and  the  problem  is  one-dimensional.  Patera  (1984)  observes  similar  behaviour 
for  spectral  element  approximations  to  second-order  problems.  The  power  and  usefulness 
of  a  multi-domain  approach  for  pseudospectral  methods  will  be  demonstrated  for  problems 
defined  in  nonrectangular  geometries  in  2-D. 

7.2  2-D  Problems 

Numerical  solutions  to  the  biharmonic  equation  are  obtained  using  the  pseudospectral  method 
when  the  exact  solution  is  given  by 

(a)  0(x,y)  =  (1  -  x’)^(l  -  y*)’sin(Ty), 

(b)  0(x,y)  =  (1  -  x^)’(l  -  y’)^3m(xx)sm(»-y), 

(c)  ^'(x.y)  =  s«r»(2xx)sm(2»y). 

In  examples  (a)  and  (b)  the  boundary  conditions  are  homogeneous  whereas  for  (c)  the 
Neumann  boundary  condition  is  inhomogeneous.  The  mixed  second  order  derivative 
is  zero  at  the  four  corners  of  fl  for  these  three  model  problems.  The  biharmonic  equation 
is  collocated  at  the  Cartesian  product  of  generalized  Legendre  nodes.  The  weighted  and 
infinity  norms  of  the  errors  are  shown  in  Table  9  for  problems  (b)  and  (c).  We  see  that  a 
numerical  solution  correct  to  machine  accuracy  is  obtained  on  a  grid  as  coarse  as  21  x  21. 
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In  the  case  of  domain  decomposition  we  consider  a  partition  of  fl  into  two  rectangular 
subdomains  (—1,0)  x  [—1,1]  and  (0, 1)  x  (—1,1)  with  common  interface  x  =  0.  Here  we  solve 
the  model  problem  (a)  using  the  collocation  scheme  (84)-(89).  Both  the  mono-domain  and 
two-domain  pseudospectral  approximations  converge  exponentially  as  expected.  Again  the 
two-domain  approximation  converges  slower  than  the  mono-domain  approximation  for  the 
same  total  number  of  collocation  points.  This  phenomenon  was  observed  for  1-D  problems 
too.  In  Figure  1  we  give  the  contours  of  the  approximation  to  the  solution  of  problem  (a) 
obtained  using  domain  decompjosition  with  N  =  20.  This  figure  is  included  to  show  the 
smoothness  of  the  contours  across  the  interface  x  =  0. 

7.3  2-D  Problems  in  Nonrectangular  Domains 

We  extend  the  ideas  developed  in  this  paper  to  the  solution  of  the  Stokes  problem  for  the 
flow  through  an  L-shaped  channel.  The  flow  geometry  in  this  example  is  nonrectangular 
for  which  a  standard  single  domain  pseudospectral  approximation  is  not  applicable.  The 
ability  of  pseudospectral  methods  to  solve  problems  in  this  kind  of  geometry  justifies  the 
development  of  the  theory  of  the  multi-domain  formulation  considered  earlier.  The  flow 
domain  is  divided  into  three  rectangular  subdomains  as  shown  in  Fig.  2.  The  stream 
function  within  each  subdomain  is  approximated  by  a  pseudospectral  representation  which 
interpolates  values  of  the  stream  function  at  interior  collocation  points  and  values  of  the 
stream  function  and  its  normal  derivative  on  the  boundaries  and  subdomain  interfaces.  These 
representations  are  C*  continuous  across  the  subdomain  interfaces.  The  unknowns  in  the 
pseudospectral  approximations  are  determined  from  the  collocation  scheme  derived  from 
the  discrete  variational  formulation.  This  scheme  results  in  continuous  approximations 
asymptotically. 

If  approximations  of  degree  N  aae  used  in  each  direction  in  each  subdomain  then  the 
collocation  equations  yield  a  system  of  {3N  —  S){N  —  3)  equations  for  the  {3N  —  5)(A^  —  3) 
unknowns.  A  total  of  2(  Af  —  3)  of  these  unknowns  represent  the  values  of  the  normal  deriva¬ 
tives  of  V’  at  the  interior  nodes  along  the  interfaces  between  subdomains  Qi  and  CI2 
between  subdomains  Qj  and  Qa-  The  remaining  unknown  values  are  the  nodal  values  of  0 
at  the  interior  and  interface  points  of  subregions  Hi,  and  Ha.  The  collocation  equations 
give  rise  to  a  linear  algebraic  system  Au  =  b.  The  vector  u  contains  the  nodal  values  of  rl> 
and  also  the  normal  derivative  of  0  at  the  interface  nodes.  The  block  tridiagonal  structure 
of  the  matrix  A  for  the  L-shaped  domain  is  shown  in  Fig.  3.  This  system  is  solved  using  a 
Crout  factorization  subroutine  from  the  NAG  Library  (1988).  A  more  efficient  direct  solution 
technique  which  takes  account  of  the  inherent  matrix  structure  is  the  almost  block  diagonal 
solver  of  Brankin  and  Gladwell  (1990)  which  has  been  used  in  spectral  calculations  by  Kara- 
georghis  and  Phillips  (1990,1991).  However,  this  subroutine  has  not  yet  been  incorporated 
into  the  present  algorithm. 

The  entry  and  exit  lengths,  a  and  6  respectively,  are  chosen  to  be  long  enough  to  obtain 
fully  developed  flow.  In  Figs.  4  and  5  we  show  the  contours  of  the  stream  function  for 
iV  =  14,  6  =  7,  c  =  1  with  a  =  — 3  and  o  =  —5,  respectively.  A  small  weak  vortex  is 
observed  in  the  salient  corner.  Fully  developed  flow  is  reached  within  a  channel  width  of  the 
reentrant  corner. 
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8  Conclusions 


Pseudospectral  approximations  to  the  solution  of  fourth-order  elliptic  partial  differential 
equations  are  constructed  using  a  collocation  procedure  based  on  the  nodes  of  generalized 
Gaussian  quadrature  rules.  Analytic  expressions  for  the  weights  appearing  in  these  quadra¬ 
ture  rules  are  derived  and  their  forms  for  the  generalized  Legendre  and  Chebyshev  rules  are 
given.  The  equivalence  between  a  discrete  variational  form  of  the  differential  problem  with 
suitably  defined  inner  products  and  a  collocation  scheme  is  demonstrated  when  the  collo¬ 
cation  points  are  chosen  to  be  the  zeros  of  certain  ultraspherical  polynomials.  The  natural 
choice  of  collocation  points  for  fourth-order  problems  differs  from  the  choice  for  second- 
order  problems,  viz.  the  Gauss- Lobatto  points.  The  usual  convergence  properties  of  spectral 
approximations  are  observed. 

A  domain  decomposition  procedure  based  on  the  generalized  Gauss- Legendre  nodes  is 
considered.  Pseudospectral  approximations  which  are  automatically  C’-  continuous  at  the 
subinterval  interfaces  are  used  to  represent  the  solution.  An  examination  of  the  correspond¬ 
ing  discrete  variational  problem  results  in  an  equivalent  collocation  method.  The  resulting 
approximation  is  shown  to  be  C^-  continuous  at  the  interfaces  asymptotically,  i.e.  as  the 
order  of  the  approximations  is  increased  in  each  subinterval.  The  scheme  is  amaJyzed  and  an 
error  estimate  is  derived  for  the  domain  decomposed  problem. 

For  fourth-order  problems  in  two  dimensions  we  propose  using  a  tensor  product  of  the 
one-dimensional  basis  functions  to  represent  the  solution.  The  equivalence  between  the  collo¬ 
cation  method  defined  by  collocating  the  differential  equation  on  a  grid  formed  by  the  tensor 
product  of  the  one-dimensional  collocation  points  and  a  discrete  variational  formulation  of 
the  problem  is  described  as  well  as  the  corresponding  domain  decomposition  problem.  It 
is  intended  to  apply  this  collocation  method  to  the  solution  of  the  Navier-Stokes  equations 
in  rectangularly  decomposable  domains  using  a  stream  function  formulation  even  though  a 
simple  variational  principle  does  not  exist  for  these  equations. 

An  application  of  this  methodology  to  a  biharmonic  problem  in  a  nonrectangular  geom¬ 
etry  is  described.  A  single  domain  approach  is  not  feasible  for  this  cl2iss  of  problems  unless 
one  first  transformed  the  original  irregular  domain  to  a  simpler  rectangular  one.  However, 
this  would  be  cumbersome  if  it  could  be  done  at  all  since  a  transformation  would  need  to  be 
found  for  each  new  geometry. 


1 


22 


TABLE  1 


Extreme  eigenvalue  of  A 


N 

■^min 

wsmstsM 

8 

0.3128+2 

0.4768+4 

0.2842-3 

12 

0.3128+2 

0.1091+6 

0.2537-3 

16 

0.3128+2 

0.1111+7 

0.2587-3 

20 

0.3128+2 

0.6788+7 

0.2652-3 

24 

0.3128+2 

0.2979+8 

0.2706-3 

28 

0.3128+2 

0.1039+9 

0.2750-3 

32 

0.3128+2 

0.3065+9 

0.2788-3 

TABLE  2 

Extreme  eigenvalue  of  A 


N 

^fnin 

^max 

8 

1.000 

3.312 

12 

1.000 

4.180 

16 

1.000. 

4.635 

20 

1.000 

4.915 

24 

1.000 

5.104 

28 

1.000 

5.241 

32 

1.000 

5.344 

TABLE  3 

Quadrature  error  in  the  approximation  of  /*,  u>(i)e*c<M(xi)dx 
for  different  weight  functions 


N 

KIOBI 

rnsjamsam 

5 

0.497  -2 

0.689  -2 

9 

0.767  -10 

0.122  -8 

17 

0.300  -15 

0.710  -14 

TABLE  4 

Quadrature  error  in  the  approximation  of  w{x)xe*dx 
for  different  weight  functions 


N 

mnmm 

msmasmM 

5 

0.579  -5 

0.843  -5 

9 

0.800-14 

0.640  -14 

17 

0.300  -15 

0.360  -14 

23 


TABLE  5 

Errors  in  the  numerical  solution  of  the  model  problem  (1) 
with  exact  solution  given  by  u{x)  =  (1  —  x^)^sin(7rx) 


N 

II  = 

HH 

IB89 

■BH 

8 

1.387-2 

1.515-2 

12 

2.941-6 

2.954-6 

16 

3.077-10 

3.041-10 

6.355-9 

24 

1.177-13 

1.329-13 

1.804-13 

2.018-13 

TABLE  6 

Errors  in  the  numerical  solution  of  the  fourth-order  problem 
with  u(±l)  =  1,  du/dx(±l)  =  2ir 
and  exact  solution  «(i)  =  1-1-  sm(27rx) 


N 

iHiig 

jRBiJI 

IIHIB 

luaui 

IBM 

■BlHi 

8 

0.266 

0.494 

12 

1.883-4 

1.845-4 

16 

1.181-7 

24 

5.795-13 

idwilpl 

4.534-13 

5.190-13 

TABLE  7 

Errors  in  the  numerical  solution  of  the  model  problem  (1)  with  exact  solution  given  by 
u(i)  =  (I  —  x^)^3in(xx)  using  domain  decomposition 


/V 

malm 

IB9 

Dsnsi 

■BIM 

HBHia 

■BHi 

6 

2.633-2 

2.024-2 

3.294-2 

2.747-2 

8 

1.096-3 

1.890-3 

1.377-3 

12 

3.494-7 

2.642-7 

4.075-7 

3.073-7 

16 

1.293-11 

1.243-11 

1.254-11 

9.523-11 

24 


TABLE  8 


Errors  in  the  numerical  solution  of  the  fourth-order  problem 
with  «(±1)  =  1,  du/dx(±l)  =  2ir  and  exact  solution  u(x)  =  1  -b  stn(2irx)  using  domain 

decomposition 


N 

non 

■iDia 

Euaia 

■BIH 

HIM 

6 

0.326 

0.476 

8 

3.580-2 

2.632-2 

12 

3.757-5 

2.839-5 

16 

7.491-9 

5.677-9 

7.%l-9 

TABLE  9 

Errors  in  the  numerical  solution  of  the  biharmonic  problem  (63)-(65) 
with  exact  solutions  given  by  (b)  and  (c) 


yv 

II  e  II3.W 

■na 

P.347 

1.467 

2.483 

8 

5.145-3 

Mn^RiVM 

6.831-2 

12 

4.667-6 

1.879-4 

4.925-4 

16 

5.631-10 

7.83e-9 

3.856-8 

7.437-8 

20 

2.935-13 

4.291-12 

9.648-11 

1.502-10 

25 


Figure  1.  Contour  plots  of  problem  (n)  when  N  =20,  using  domain 

decompodtion  and  the  generalised  Legendre  psendospectral  method. 


Figare  5.  Contour*  of  f®*"  ®  =  “®*  6  =  7,  c  =  1 
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A  Properties  of  Ultraspherical  Polynomials 

The  ultraspherical  or  Gegeobauer  Polynomials  are  the  solutions  of  the  differential  equation 


[u>a+i(x)G^‘'>'(x)]'  +  n(n  +  2A  +  l)u;A(x)Gj,^Hx)  =  0, 
that  are  bounded  at  x  =  ±  1,  where 

u;a(x)  =  (1  -  x^)^,  A  >  -1. 

They  are  orthogonal  with  respect  to  wx(x)  over  the  interval  [-1,1]  : 

tOA(x)GiJHx)Gl^>{x)dx  =  'tmSm.n  , 


where 


7»n  — 


+  A  + 1)]^ 


(2m  -I-  2A  +  l)m!r{m  -I-  2A  -|- 1)  ’ 
and  r  is  the  gamma  function.  At  x  =  ±1  ,  (?^^l(x)  satisfies  the  condition 

G<-^^(±1)  =  (±1)"^^”— — 

"  ^  ^  ^  ^  n!r(A-|-l) 

The  ultraspherical  polynomials  may  be  generated  using  the  recurrence  relation 

(n+l)(n  +  2A  +  l)C!<4>. 

=  (2r.  +  2A  +  l)(n  +  A  +  l)iG!,»>  -  (n  +  A)(n  +  A  +  1)G?.>„ 
G<o*’(i)=l,  Gi‘'(i)  =  (A+l)i. 

The  leading  coefficient,  An,  of  is  given  by 

1  r(2n  +  2A-|-l) 

Ay^  — 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


2"  n!r(n  d- 2A -i- 1)  ■ 

We  have  the  following  integrals  involving  ultraspherical  polynomials  (Erdelyi  (1954),  p.284) 


»!r(A  -  (r)r(A 

where  A,  <r  >  — 1. 

The  ultraspherical  polynomials  satisfy  the  recursion  relation 

(1  -  x”)G<;)'(x)  =  -nxG<,*>(x)  +  (o  +  A)G<,'J,(x). 


(10) 
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